* Get PREDICTED wind capacity around a county, as an IV for actual wind capacity around a location: macro; 

* define libraries;
libname hannah  '/projects/users/########';
libname ben '/projects/users/########';
libname hantmp  '/projects/users/########/tempfiles';

proc import out = work.cnty_level_ivs2
  datafile = "/projects/users/########/cnty_level_ivs2.dta"
  dbms = dta
  replace;
run;

data work.cnty_level_ivs2;
  set work.cnty_level_ivs2 
  (keep = ctyfips year zhat_cnty cond_exp_ycap cond_exp_ytnum meanspd);
run;

proc sql;
  create table work.cnty_level_ivs2 as select
  a.*, b.x, b.y
  from work.cnty_level_ivs2 as a
    left join hannah.countylevelwindwrkrs as b on a.ctyfips = b.ctyfips and a.year = b.year
  order by a.ctyfips, a.year;
quit;

** Now loop over the circles to get PREDICTED wind capacity at each distance from each county.
*define macro "cnty_iv" which does X-mile county buffer for all years, counties;
%macro cnty_iv;

  *get centers of all counties;
  data work.mile20centers;
    merge hannah.USmile20circlesc hannah.centersc;
    by _ID_;
  run;
  data work.mile40centers;
    merge hannah.USmile40circlesc hannah.centersc;
    by _ID_;
  run;
  data work.mile60centers;
    merge hannah.USmile60circlesc hannah.centersc;
    by _ID_;
  run;
  data work.mile80centers;
    merge hannah.USmile80circlesc hannah.centersc;
    by _ID_;
  run;
  data work.mile100centers;
    merge hannah.USmile100circlesc hannah.centersc;
    by _ID_;
  run;

  %do j = 1 %to 1417;
    *get individual county buffers; 
    data work.usmile20_&j;
      set hannah.USmile20circlesc (keep = X Y SEGMENT _ID_);
      if _ID_ = &j;
    run;

    data work.usmile40_&j;
      set hannah.USmile40circlesc (keep = X Y SEGMENT _ID_);
      if _ID_ = &j;
    run;

    data work.usmile60_&j;
      set hannah.USmile60circlesc (keep = X Y SEGMENT _ID_);
      if _ID_ = &j;
    run;

    data work.usmile80_&j;
      set hannah.USmile80circlesc (keep = X Y SEGMENT _ID_);
      if _ID_ = &j;
    run;

    data work.usmile100_&j;
      set hannah.USmile100circlesc (keep = X Y SEGMENT _ID_);
      if _ID_ = &j;
    run;

    *get individual county centers;
    data work.uscntycenter_&j;
      set hannah.centersc (keep = _ID_ meanx meany);
      if _ID_ = &j;
    run;
      
      *GET PLANTS WITHIN THE BUFFER OF &j;
      *Create temp file of all counties near county j in any year;
      proc sql;
	create table work.cnties_&j as
	select 
	  A.*
	from work.cnty_level_ivs2 as A
	  left join work.uscntycenter_&j as E on ctyfips
	where A.y between E.meany + 3 and E.meany - 3
	  and A.x between E.meanx + 3 and E.meanx - 3
	  and GEODIST(A.y, A.x, E.meany, E.meanx, 'M' ) <= 103
	  and GEODIST(A.y, A.x, E.meany, E.meanx, 'M' )  ^=.
	order by A.ctyfips, A.year;
      quit;
      
      *Flag observations inside X miles, 20, 40, 60, 80, 100;
      PROC GINSIDE 
	DATA= work.cnties_&j
	MAP = work.usmile20_&j
	OUT = hannah.site_Xm_&j;
	ID _ID_;
      RUN;
      data hannah.site_Xm_&j;
	set hannah.site_Xm_&j;
	IF CMISS(_ID_) THEN _ID_ = 0;
	rename _ID_ = id20;
      run;
      PROC GINSIDE 
	DATA= hannah.site_Xm_&j
	MAP = work.usmile40_&j
	OUT = hannah.site_Xm_&j;
	ID _ID_;
      RUN;
      data hannah.site_Xm_&j;
	set hannah.site_Xm_&j;
	IF CMISS(_ID_) THEN _ID_ = 0;
	rename _ID_ = id40;
      run;
      PROC GINSIDE 
	DATA= hannah.site_Xm_&j
	MAP = work.usmile60_&j
	OUT = hannah.site_Xm_&j;
	ID _ID_;
      RUN;
      data hannah.site_Xm_&j;
	set hannah.site_Xm_&j;
	IF CMISS(_ID_) THEN _ID_ = 0;
	rename _ID_ = id60;
      run;
      PROC GINSIDE 
	DATA= hannah.site_Xm_&j
	MAP = work.usmile80_&j
	OUT = hannah.site_Xm_&j;
	ID _ID_;
      RUN;
      data hannah.site_Xm_&j;
	set hannah.site_Xm_&j;
	IF CMISS(_ID_) THEN _ID_ = 0;
	rename _ID_ = id80;
      run;
      PROC GINSIDE 
	DATA= hannah.site_Xm_&j
	MAP = work.usmile100_&j
	OUT = hannah.site_Xm_&j;
	ID _ID_;
      RUN;
      * keep id100 named _ID_;

      * Drop observations outside the buffer;
      DATA hannah.site_Xm_&j;
	SET hannah.site_Xm_&j;
	IF CMISS(_ID_) THEN DELETE;
	plntid = &j;
      RUN;

      data hannah.site_Xm_&j;
	merge hannah.site_Xm_&j (IN = thisj) hannah.centersc;
	by _ID_;
	if thisj = 1;
	drop _TYPE_ _FREQ_;
      run;
      * Collapse to the county-year level; 
      proc sql;
	CREATE TABLE work.cnties_sum_&j AS
	SELECT _ID_, year,
	SUM(zhat_cnty) as zcnty100,
	SUM(case when id80~=0 then zhat_cnty else 0 end) as zcnty80,
	SUM(case when id60~=0 then zhat_cnty else 0 end) as zcnty60,
	SUM(case when id40~=0 then zhat_cnty else 0 end) as zcnty40,
	SUM(case when id20~=0 then zhat_cnty else 0 end) as zcnty20,
	SUM(cond_exp_ycap) as yhatcap100,
	SUM(case when id80~=0 then cond_exp_ycap else 0 end) as yhatcap80,
	SUM(case when id60~=0 then cond_exp_ycap else 0 end) as yhatcap60,
	SUM(case when id40~=0 then cond_exp_ycap else 0 end) as yhatcap40,
	SUM(case when id20~=0 then cond_exp_ycap else 0 end) as yhatcap20,
	SUM(cond_exp_ytnum) as yhattnum100,
	SUM(case when id80~=0 then cond_exp_ytnum else 0 end) as yhattnum80,
	SUM(case when id60~=0 then cond_exp_ytnum else 0 end) as yhattnum60,
	SUM(case when id40~=0 then cond_exp_ytnum else 0 end) as yhattnum40,
	SUM(case when id20~=0 then cond_exp_ytnum else 0 end) as yhattnum20,
	MEAN(meanspd) as mwspeed100,
	MEAN(case when id80~=0 then meanspd else . end) as mwspeed80,
	MEAN(case when id60~=0 then meanspd else . end) as mwspeed60,
	MEAN(case when id40~=0 then meanspd else . end) as mwspeed40,
	MEAN(case when id20~=0 then meanspd else . end) as mwspeed20,
	MAX(ctyfips) as ctyfips,
	MAX(x) as mx,
	MAX(y) as my
	FROM hannah.site_Xm_&j
	GROUP BY _ID_, year
	ORDER BY _ID_, year;
      quit;
      *get individual location ctyfips; 
      data work.usmileXvars_&j;
	set work.mile20centers;
	if _ID_ = &j;
      run;
      proc sql;
	create table work.cnty_&j as
	select _ID_,
	MAX(ctyfips) as ctyfips
	FROM work.usmileXvars_&j
	GROUP BY _ID_
	ORDER BY _ID_;
      quit;
      data work.cnties_sum_&j;
	merge work.cnties_sum_&j work.cnty_&j;
	by _ID_;
      run;
	
      PROC APPEND base = hannah.cnty_dist_nearIV
	data = work.cnties_sum_&j;
      RUN;
      PROC DATASETS library = work nolist;
	delete cnties_sum_&j cnty_&j usmileXvars_&j cnties_&j;
      RUN; QUIT;
      PROC DATASETS library = hannah nolist;
	delete site_Xm_&j;
      RUN; QUIT;

    PROC DATASETS library = work nolist;
      delete uscntycenter_&j usmile20_&j usmile40_&j usmile60_&j usmile80_&j usmile100_&j;
    RUN; QUIT;
  %end;
  PROC DATASETS library = work nolist;
    delete mile20centers mile40centers mile60centers mile80centers mile100centers;
  RUN; QUIT;
%mend;

/*
proc means noprint data=hannah.USmile100circlesc;
  var X Y;
  by _ID_;
  output out = hannah.centersc mean(X Y) = meanx meany;
run;


PROC DATASETS library = hannah nolist;
  delete cnty_dist_nearIV;
RUN; QUIT;
*/

%cnty_iv; 

/*
data hannah.cnty_shp_ids;
  set hannah.usmile40circlesc (keep = _ID_ ctyfips);
run;
proc sort data = hannah.cnty_shp_ids nodupkey;
  by _ID_ ctyfips;
run;
data hannah.cnty_dist_nearIV;
  set hannah.cnty_dist_nearIV (drop = ctyfips);
run;
proc sql;
  create table hannah.cnty_dist_nearIV as select
  a.*, b.ctyfips
  from hannah.cnty_dist_nearIV as a
    left join hannah.cnty_shp_ids as b
    on a._ID_ = b._ID_
  order by a._ID_, a.year;
quit;
proc sort data= hannah.cnty_dist_nearIV;
  by ctyfips year;
run;
*/
